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Abstract. 

A new Chemo - Dynamical Smoothed Particle Hydrodynamic (CD - 
SPH) code is presented. The disk galaxy is described as a multi - frag- 
mented gas and star system, embedded in a cold dark matter halo. The 
star formation (SF) process, SNII, SNIa and PN events as well as chem- 
ical enrichment of gas have been considered within the framework of the 
standard SPH model. Using this model we describe the dynamical and 
chemical evolution of triaxial disk - like galaxies. It is found that such ap- 
proach provides a realistic description of the process of formation, chemical 
and dynamical evolution of disk galaxies over a cosmological timescale. 



1. Introduction 

The dynamical and chemical evolution of galaxies is one of the most inter- 
esting and complex problems. Naturally, galaxy formation is tightly con- 
nected with the process of large - scale structure formation in the Universe. 
The main role in the scenario of large - scale structure formation seems to 
be played by the dark matter. It is believed that the Universe was seeded 
at some early epoch with small density fluctuations of dark non - baryonic 
matter, and the evolving distribution of these dark halos provides the arena 
for galaxy formation. Galaxy formation itself involves collapse of baryons 
within potential wells of dark halos (White & Rees, 1978). The properties 
of forming galaxies depend on the amount of baryonic matter, that can be 
accumulated in such halos, and the efficiency of star formation. The obser- 
vational support for this scenario of galaxy formation comes from the recent 
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COBE detection of fluctuations in the microwave background e.g. (Bennett 
et ai, 1993). 

The investigation of the process of galaxy formation is a highly complex 
subject requiring many different approaches. The formation of self - gravi- 
tating inhomogeneities of protogalactic size, the ratio of baryonic and non 
- baryonic matter (Bardeen et ai, 1986; White & Silk, 1979; Peebles, 1993; 
Dar, 1995) the origin of the protogalaxy's initial angular momentum (Voglis 
&: Hiotelis, 1989; Zurek et al., 1988; Eisenstein & Loeb, 1995; Steinmetz & 
Bartelmann, 1995), the protogalaxy's collapse and its subsequent evolution 
are usually considered as separated problems. The recent advances in com- 
puter technology and numerical methods have made possible the detailed 
modeling of baryon matter dynamics in the universe dominated by colli- 
sionless dark matter and, therefore, the detailed gravitational and hydro- 
dynamical description of the formation and evolution of galaxies. The most 
sophisticated models additionally include radiative processes, star forma- 
tion and supernova feedback e.g. (Katz, 1992; Steinmetz & Muller, 1994; 
PriedU & Benz, 1995). 

However, results of numerical simulations are essentially affected by 
the star formation algorithm incorporated into modeling techniques. The 
star formation and accompanied processes are still not well understood 
on either small or large spatial scales. The star formation algorithm by 
which the gas material is converted into stars can be based only on simple 
theoretical assumptions or on results of empirical observations of nearby 
galaxies. The other most important effect of star formation on the global 
evolution of a galaxy is caused by a large amount of energy released in 
supernova explosions and stellar winds. 

Among numerous methods developed for the modeling of complex three 
dimensional hydrodynamic phenomena Smoothed Particle Hydrodynamics 
(SPH) is one of the most popular (Monaghan, 1992). Its Lagrangian nature 
allows to be combined easily with fast N - body algorithms, making it very 
suitable for simultaneous description of complex dynamics of gas - stellar 
system (Fricdli & Bcnz, 1995). As an example of such combination, the 
TREE - SPH code (Hernquist & Katz, 1989; Navarro k White, 1993) can 
be mentioned, which was successfully applied to the detailed modelling 
of disk galaxy mergers (Mihos Sz Hernquist, 1996) and galaxy formation 
and evolution (Katz, 1992). The second good example is a GRAPE - SPH 
code (Steinmetz &: Muller, 1994; Steinmetz & Muller, 1995) which was 
successfully used to model of the evolution of disk galaxy structure and 
kinematics. 

In recent few years, we had an high number of excellent papers concern- 
ing the complex SPH modeling of galaxy formation and evolution (Raiteri 
et ai, 1996; Carraro et al., 1997). In our code, we propose new "energetic" 
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criteria for SF, and make a more realistic account of returned chemically 
enriched gas fraction via SNII, SNIa and PN events. 

The simplicity and numerical efficiency of the SPH method were the 
main reasons why we chose this technique for the modeling of the evolution 
of complex, triaxial protogalactic systems. We used our own modification 
of the hybrid N - body/SPH method (Bcrczik & Kravchuk, 1996; Berczik, 
1998), which we call the chemo - dynamical SPH (CD - SPH) code. 

The "stars" are included into the standard SPH algorithm as the N 
- body collisionless system of particles, which can interact with the gas 
component only through the gravitation (Katz, 1992). The star formation 
process and supernova explosions are included into the scheme in the man- 
ner proposed by Raiteri et al. (1996), but with our own modifications. 

2. The CD - SPH code 

2.1. THE SPH CODE 

Continuous hydrodynamic fields in SPH are described by the interpolation 

functions constructed from the known values of these functions at randomly 
positioned particles. In this case the mean value of a physical field /(r) at 
the point r can be written as (Monaghan, 1992): 

</(r)>= I f{r').W{r-r';h)dr', (1) 

where W{r — r'; h) is the kernel function and h is the softening constant. 

The function W{r — r'; h) is strongly peaked at | r — r' |= 0, and we can 
consider it without any loss of generality as belonging to the class of even 
functions. In this case it is not difficult to demonstrate (Hernquist & Katz, 
1989) that the average value < /(r) > represents the real function /(r), 
with an error not higher than 0{h^). Consider a fluid with the density /^(r), 
then we can rewrite (1) in the form: 

< /(r) >= / ^ ■ W^(r - r'; h) ■ p{v')dv' . (2) 

Let us imagine that /(r) and p(r) are known only at A'" discrete points 
Tj. Then, equation (2) gives: 

N r 

fi = T.^r^■W^r (3) 

Here = /(r^), pi = p(rj), Wij = W{ri - Tf,h) and m.^ is the mass 
of particle i. Using (3), we approximate the hydrodynamic field by an an- 
alytical function which is differentiable as many times as the kernel Wy. 
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Following Monaghan & Lattanzio (1985) we use for the kernel function Wij 
the spline expression in the form: 



1 - + f ufj-, if < Uij < 1, 

1{2-Uijf, iil<Uij<2, (4) 

0, otherwise . 



Here uij = Vij/h. 

The local resolution in SPH is defined by the chosen smoothing length h. 
If it is the same for all points, the hydrodynamic field will be evidently ap- 
proximated more smoothly in the regions where the points lie more closely, 
i.e. where the density is higher. 

To achieve the same level of accuracy for all points in the fluid it is 
necessary to use a spatially variable smoothing length. In this case each 
particle has its individual value of h. Following Hernquist &; Katz (1989), 
we write instead of (1): 

< /(r) >= j f{v') ■ ^ ■ [W{5v- h) + W{5v- h')]dv' . (5) 

Here (5r = r — r' , /i = /i(r) and h! = h{r'). For the density p{r) it gives 

AT ^ 

Other hydrodynamic functions are written in the same manner. 

In our calculations the values of hi were determined from the condition 
that the number of particles Nb in the neighborhood of each particle within 
the 2 • hi remains constant (Mihos & Hernquist, 1996). The value of Nb is 
chosen such that a certain fraction from the total number of "gas" particles 
N affects the local flow characteristics (Hiotclis &; Voglis, 1991). If the 
defined hi become smaller than the minimal smoothing length hmin, we set 
the value hi = hmin- For "star" particles (with Plummer density profiles) 
we use, accordingly, the fixed gravitational smoothing lengths hgtar- 

If the density is computed according to equation (6), then the continuity 
equation is satisfied automatically. Equations of motion for particle i are 

^ = + _ v,$, - Vi^-\ (8) 

(It Pi 
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where Pi is the pressure, $i is the self gravitational potential, <&?^* is a 
gravitational potential of possible external halo, aj"** is an artificial viscosity 
term. Using equations (5) and (6) we can write: 



v.P. ^ ^ ^mw{nr,h,) + Winy,h,)]. (9) 



Pi 



=1 



The artificial viscosity term a^*^ is introduced so as to describe the 
flows with shock waves. In present calculations the viscous acceleration is 
introduced by replacing [Pi j + Pj / p^) in equation (9) by {Pi / pf + Pj / pj) ■ 
(1 + TTij). The expression for TTy has the form (Hiotelis et al, 1991): 



-a ■ Hij + P ■ 



(10) 



where a and /3 are constants, and Hij is defined by the relation: 



f^ij 



hij/cij ■ {^ijVij)/{r1j + n^hjj), if (vyTy) < 0, 
0, otherwise. 



(11) 



Here Cij = {ci + Cj)/2 is the sound speed averaged for points i and 
J, hij = {hi + hj)/2, Yij = {vi — Tj) and \ij = (vj — Vj) is the relative 
velocity vector for the points i and j. The term • hf^ is inserted to prevent 
divergences when Vij = and it should be small enough. The constant n 
was set equal to 0.1. For the constant a and (3 in (10) the values a = 1 and 
j3 = 2 give good results in a wide range of the Mach numbers (Monaghan, 
1992). 

By using equation (4), one gets for the gravitational acceleration (Hiotelis 
& Voglis, 1991): 



1 



N 



(12) 



where: 



9k 



4, ,3 _ 6 5 , 1,,6 



3'"fc 5"*; 



< ■Ufc < 1, 
otherwise. 



(13) 



and Ujfc = rij/hk- 
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When isothermal flows are considered, the system of equations is closed 
by adding the equation of state: 

P^ = P^■cl (14) 

where Cj is the isothermal speed of sound. 

In the particle representation for complex hydrodynamic flows the en- 
ergy equation has the form: 



X Vi[P^(ri,-; hi) + W{nj- hj)\wij + (15) 

Pi 

Here ui is the specific internal energy of particle i. The term (Fj — 

Ki)/ Pi accounts for nonadiabatic processes not associated with the artificial 
viscosity (in our calculations Fj = 0). We present the radiative cooling in 
the form: 

ki = Ki{ui,pi)=k*{Ti)-nj, (16) 

where rii is the hydrogen number density and Tj is the temperature. 
To follow its subsequent thermal behaviour in numerical simulations we 
use an analytical approximation of the standard cooling function A* (T) for 
an optically thin primordial plasma in ionization equilibrium (Dalgarno &; 
McCray, 1972; Katz k, Gunn, 1991). Its absolute cutoff temperature is set 
equal to 10^ K. 

If log r < 6.2: 

A*(r) 



else 



10 



-28 



+ 10-1.7-0.2(6.2-logTr] j/g^ 



(17) 



A*(T) = 10"^^ • [io-0-i-i-88(5-23-iogr)4 _^ 

+10-1-^] J/s, (18) 

The equation of state must be added to close the system. 



Pi = Pi-{l ■ Ui, 



(19) 
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where 7 = 5/3 is the adiabatic index. 



2.2. TIME INTEGRATION 

To solve the system of equations (7), (8) and (15) we use the leapfrog 
integrator (Hernquist & Katz, 1989). This system of equations has the 
form: 



' dxi/dt = Vi, 
dvi/dt = SLi{h, p, T, c,P,r,v), 



(20) 



, dui/dt = Ui{h,p,T,c,P,r,v). 

The time step Sti for each particle depends on the particle's acceleration 
aj and velocity Vj, as well as on viscous forces. To define 5ti we use the 
relation (Hiotelis & Voglis, 1991): 



Sti = Cn min[< 



hi 



hi hi 



I ' J ' 

Vi Si 



(21) 



where Si = Ci ■ {1 + a + 0.68 • P ■ maxj \ fiij |) and Cn is the Courant's 
number. We adopt Cn = 0.1. The minimum time step value is set equal to: 



At = min (Sti), 



(22) 



for all particles. The integrator has a second order accuracy in the time 
step At. As the first step for particle i we define: 



' r"+^ =r2 + 0.5 • At ■ wf, 

< v^''+^ = + 0.5 • At • af , (23) 

n+i 

_ =uf + 0.5-At-uf. 

After that: 



' r^+i = rf + At- v"+^ + 0{At% 

< v^+i = vf + At-a"+^ +0(Ai3), (24) 
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We have carried out (Berczik & Kolesnik, 1993) a large series of test 
calculations to check that the code is correct, the conservation laws are 
obeyed and the hydrodynamic fields are represented adequately. These tests 
have shown very good results. 

2.3. THE STAR FORMATION ALGORITHM 

It is well known that star formation (SF) regions are tightly connected to gi- 
ant molecular complexes, especially to regions, which reach some threshold 
for dynamical instabilities. The overall picture of star formation seems to 
be understood, but the detailed physics of star formation and accompany- 
ing processes on either small or large scales remains sketchy (Larson, 1969; 
Silk, 1987; Leitherer et al, 1992). 

For describing of the process of conversion of gas material into stars we 
modify the standard SPH star formation algorithm (Katz, 1992; Navarro 
& White, 1993), taking into account the presence of chaotic motions in the 
gaseous environment and the time lag between the initial development of 
suitable conditions for star formation and star formation itself (Berczik &; 
Kravchuk, 1996; Berczik, 1998). The first reasonable requirement incorpo- 
rated into this algorithm allows selecting "gas" particles that are potentially 
eligible to form stars. It states that in the separate "gas" particle the SF 
can start if the absolute value of the "gas" particle's gravitational energy 
exceeds the sum of its thermal energy and energy of chaotic motions: 

\Ef \> Ef' + Ef. (25) 

Gravitational and thermal energies and energy of chaotic motions for 
the "gas" particle i in model simulation are defined as: 

' Er = -I.G.mVhi, 

< £;f = |-mi-c2, (26) 
_ Eff' = ^-mi- Avl 

Where Cj = ^/W^Tifjl is the isothermal sound speed of particle i. We 
set /X = 1.3 and define the chaotic or "turbulent" square velocities near 
particle i as: 

Nb Nb 
Avf = J^mj- {wj - Ve)V ^ mj, (27) 

where: 
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Nb Nb 
Vc = • Vj/^mj. (28) 

i=i i=i 

For practical reasons, it is useful to define a critical temperature for SF 
onset in particle i as: 

Tr' = ^■il■'^■G■p^■hi-i^vl). (29) 

Then if the temperature of the "gas" particle z, drops below the critical 
one, SF can proceed. 

Ti < T["K (30) 

We think that requirement (25), or in an other form (30), is only a nec- 
essary one. It seems reasonable that the chosen " gas" particle will produce 
stars only if the above condition will hold over the time interval exceeding 
its free - fall time tff = -\/3 • tt/ (32 ■ G ■ p). This condition is based on the 
well known fact that due to gravitational instability all substructures of 
collapsing system are formed on such timescale. Using it wc exclude from 
the consideration transient structures, that are destroyed by the tidal ac- 
tion of surrounding matter. This last condition we assume to be a sufficient 
one. 

We also define which "gas" particles remain cool, i.e. tcooi < ^//- We 
rewrite this condition in the manner presented in Navarro &; White (1993): 
Pi > Pcrit- Here we use the value of pcru = 0.03 cm~^. 

When the collapsing particle i is defined, we create the new "star" 
particle with mass m**"'' and update the "gas" particle mj using these 
simple equations: 



^star _ g . j^^^ 

(31) 

mj = (1 - e) ■ rrii. 

Here e, defined as the global efficiency of star formation, is the fraction of 
gas converted into stars according to the appropriate initial mass function 
(IMF). The typical values for SF efficiency in our Galaxy on the scale of 
giant molecular clouds are in the range e « 0.01 — 0.4 (Duerr et al, 1982; 
Wilking & Lada, 1983). But it is still a poorly known quantity. In numerical 
simulation it is the model parameter which has to be checked by comparison 
of numerical simulation results with available observational data. Here we 
define e as: 
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e = l-{Ef' + Ef)/ \Ef \, (32) 

with the requirement that all excess mass of the gas component in star 
forming particle, which provide the inequality | Ef^ |> Ef^ + Ef^. is trans- 
formed into the star component. In the code we set the absolute maximum 
value of the mass of such a "star" particle m^^^ = 2.5 • 10^ Mq, i.e. w 5% 
of the initial particle mass mi. 

At the moment of the birth, the position and velocities of new "star" 
particles are equal to those of parent "gas" particles. Thereafter these "star" 
particles interact with other " gas" and " star" particles only by gravitation. 
The gravitational smoothing length for these (Plummer like) particles is 
set equal to hgtar- 

2.4. THE THERMAL SNH FEED - BACK 

We try to include in our code the events of SNII, SNIa and PN into the 
complex gasdynamic picture of galaxy evolution. But, for the thermal bud- 
get of the ISM, only SNII play the main role. Following Katz (1992), we 
assume that the explosion energy is converted totally to the thermal energy. 
The stellar wind action seems not to be essential in the energy budget (Fer- 
riere, 1995). The total energy released by SNII explosions (10^"^ J per SNII) 
within a " star" particles is calculated at each time step and distributed uni- 
formly between the surrounding (i.e. Vij < hgtar) "gas" particles (Raiteri 
et al, 1996). 

2.5. THE CHEMICAL ENRICHMENT OF GAS 

Every new "star" particle in our SF scheme represents a separate, gravi- 
tationally closed star formation macro region (like a globular cluster). The 
"star" particle is characterized its own time of birth thegSF which is set 
equal to the moment of particle formation. After the formation, these par- 
ticles return the chemically enriched gas to surrounding " gas" particles due 
to SNII, SNIa and PN events. For the description of this process we use 
the approximation proposed by Raiteri et al. (1996). We consider only the 
production of ^^O and ^^Fe, and try to describe the full galactic time evolu- 
tion of these elements, from the beginning to present time (i.e. tf.yoi ~ 13.0 
Gyr). 

We use the multi - power IMF law defined by Kroupa et al. (1993). 
The distribution of stellar masses within a "star" particle of mass m**"*" is 
then: 
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r 20-9 • 



m 



m 



m 

-2.2 
-2.7 



if 0.1 < m < 0.5, 
if 0.5 < m < 1.0, 
if 1.0 < m. 



(33) 



Here m is the star mass in solar units and the normahzation constant 
is A = 0.310146. The value of A is different from Raiteri et al. (1996), 
because we set the integration limits migyj = 0.1 Mq and m^pp = 100 Mq. 

For the definition of stellar lifetimes we use the equation (Raiteri et al., 
1996): 



\ogtdead = ao{Z) - ai(Z) • logm + a2(Z) • (logm)^, 



(34) 



where i^ead is expressed in years and m is in solar unit and coefficients 
are defined as: 



' ao(Z) = 10.130 + 0.0755 • log Z - 0.0081 • (log Z)2, 

< ai(Z) = 4.4240 -I- 0.7939 -log Z + 0.1187- (log Z)2, (35) 

_ a2(Z) = 1.2620 + 0.3385 • log Z + 0.0542 • (log Zf. 

These relations are based on the calculations of the Padova group (Alongi 
et al, 1993; Bressan et al., 1993; Bertelli et al, 1994) and give a reasonable 
approximation to stellar lifetimes in the mass range from 0.6 Mq to 120 Mq 
and metalicities Z = 7- 10~^ — 0.03 (defined as a mass of all elements heavier 
than He). In our calculation, we assume following, Raiteri et al. (1996), 
that Z scales with the oxygen abundance as Z/Zq =^^0/^^00. For lower 
or higher metalicities we take the value corresponding to the extremes of 
this interval. 

We can define the number of SNH explosions inside a given "star" par- 
ticle during the time from t to t + At using a simple equation: 

AiVsNii = j ^{m)dm. (36) 

Here mdead{t) and rudeadit + At) are masses of stars that end their 
lifetimes at the beginning and at the end of the time step respectively. 
Contrary to Raiteri et al. (1996), we assume that all stars with masses 
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between 8 Mq and 100 Mq end their lives as SNII. For SNII we use the 
yields from Woosley & Weaver (1995). We use the approximation formulae 
from Raitcri et al. (1996) for defining the total ejected mass by one SNII 
- m*"*, as well as for the ejected mass of iron - and oxygen - m^- as a 
function of stellar mass (in solar unit). 



' rra*°* = 7.682 ■ 10"^ • m^ ^^^, 

< m^l = 2.802 • 10-"^ • (37) 
m?. = 4.586 ■ IQ-^ • m^ '^^i 

To take into account PN events inside the "star" particle we use the 
equation, like (36): 

AA/^PN = j ■^{m)dm. (38) 

mdead{t+^t) 

Following van den Hoek Sz Groenewegen (1997), Samland (1997) and 
Samland et al. (1997), we assume that all stars with masses between 1 Mq 
and 8 Mq end their lives as PN. We define the average ejected masses (in 
solar unit) of one PN event as (Renzini &; Voli, 1981; van den Hoek & 
Groenewegen, 1997): 



m, 



tot 



1.63, 



ml] = 0.00, 



(39) 



m 



0.00. 



The method described in Raitcri et al. (1996) and proposed in Greggio 
& Renzini (1983) and Matteuchi &: Greggio (1986) is used to account of 
SNIa . In simulations, the number of SNIa exploding inside a selected "star" 
particle during each time step is given by: 



"^dead(*) 



AA^SNIa = 



/ 



'^2{rn2)dm2. 



(40) 



mdead{t+M) 

The quantity ^2("^2) represents the initial mass function of the sec- 
ondary component and includes the distribution function of the secondary's 
mass relative to the total mass of the binary system, m^. 
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*2(m2) = m^*"'--^2- / i—f-rriB^-'^dmB, (41) 

J vriB 

where rriinf = max(2 • m2, 3 Mq) and nisup = ^-2 + 8 Mq. The value 
of normahzation constant we set, following the van den Berg & McClure 
(1994), equal to A2 = 0.16 • A. 

Then the total ejected mass (in solar unit) is (Thielemann et ai, 1986; 
Nomoto et al, 1984): 

' mlf = 1.41, 

< m^| = 0.63, (42) 
mO = 0.13. 

In summary, a new "star" particle (with metalicity Z = 10~^) with mass 
10^ Mq, during the total time of evolution tevol produces: 

AA^SNii ~ 52.5, AATpN ^ 1770.0, AA^sNia ~ 8.48. 



10 "3- 



NPN 




1 10 10 ' 10 ' 10 * 

t (Myr) 



Figure 1. The number of SNII, SNIa and PN event 
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We present in Fig. 1. the number of SNII, SNIa and PN events for this 
"star" particle. 

In Fig. 2. and Fig. 3. we present the returned masses of ''^^'Fe and ^^O, 
respectively. We can estimates the total masses (H, He, ^^Fe, ^^O) (in solar 
masses) returned to the surrounding "gas" particles, due to these processes 
as: 



^SNII 

"-SNII 

,Fe 
'SNII 



477, 


AmpN 


= 2164, 


^"^SNIa 


= 4.14, 


159, 




= 721.3, 




= 1.38, 


3.5, 


Am^lf 


= 0.000, 


^"^SNIa 


= 5.35, 


119, 


Am^N 


= 0.000, 


Am^NIa 


= 1.10. 



(43) 



2.6. THE COLD DARK MATTER HALO 

In the literature we have found some, sometimes controversial, profiles for 
the galaxies Cold Dark Matter Haloes (CDMH) (Burkert, 1995; Navarro, 
1998). For resolved structures of CDMH: phalo{r) ~ r~^-^ (Moore et al, 
1997). The structure of CDMH high-resolution N-body simulations show 
can be described by: Phalo{f) ~ f'"^ (Navarro et al., 1996; Navarro et al., 
1997). Finally, in Kravtsov et al. (1997), we find that the cores of DM 
dominated galaxies may have a central profile: phaio{f) ^ r^°'^. 

Because we concentrate our attention mainly on the problem of the 
formation and evolution (dynamical and chemical) of the disk structure in 
galaxy we using the simplest dark matter halo model. In our calculations, 
as a first approximation, it is assumed that the model galaxy halo contains 
the CDMH component with Plummer - type density profiles (Douphole &; 
Colin, 1995): 



P^aUr) = ^ • — (44) 

Therefore we can write for the external force acting on the "gas" and 
"star" particles: 

-V,$f * = -G ■ 3 • r,. (45) 

The inclusion of the dynamically evolved halo model is a next step in 
our calculations. 
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3. Results and discussion 
3.1. INITIAL CONDITIONS 

Wc test our code and demonstrate that simple assumptions lead to a rea- 
sonable galaxy model. The SPH calculations were carried out for Ngas = 
2109 "gas" particles. According to Navarro &; White (1993) and Raiteri et 
al. (1996), such number seems to be quite enough to provide a qualita- 
tively correct description of the system behaviour. Even such small num- 
ber of "gas" particles produces Ngtar = 31631 "star" particles at the 
end of the calculation. We also check the dependency of our model with 
the number of "gas" particles. We calculate the first ~ 5 Gyr evolution 
of the systems with identical initial data but with different numbers of 
Ngas = 2109; 4169; 11513 with correspondent Nb = 21; 41; 115. All three 
models show the practically identical star formation rate and chemical evo- 
lution history. 

The value of the smoothing length hi was chosen requiring that each 
"gas" particle had Nb = 21 neighbors within 2 ■ hi. Minimal hmin was set 
equal to 1 kpc. This value is defined as a largest characteristic dimensions of 
instability in galactic disk. This hmin is enough only for a crude description 
of the galaxy disk structure but is quite enough for a dimensions of SSP 
(Single Stellar Population) in chemical evolution modeling. Any dynami- 
cal processes below this scale we can't describe correctly in our SPH code. 
We use also the fixed gravitational smoothing length hstar = 1 kpc for the 
"star" particles. Our results show that such a value of Nb ~ 1% Ngas pro- 
vides qualitatively correct treatment of the system's large scale evolution. 

As the initial model (relevant for CDM - scenario) we took constant 
- density homogeneous gaseous triaxial configuration {Mgas = 10^^ Mq) 
within the dark matter halo {Mhalo = 10^^ Mq). We set A = 100 kpc, B = 
75 kpc and C = 50 kpc for semiaxes of system. Such a triaxial configurations 
are reported in cosmological simulations of the dark matter halo formation 
(Eisenstein & Loeb, 1995; Prenk et al, 1988; Warren et al, 1992). Initially, 
the centers of all particles were placed on a homogeneous grid inside this 
triaxial configuration. Wc set the smoothing parameter of CDMH: b^alo = 
25 kpc. Such values of M^alo and h^alo are typical for CDMH in disk galaxies 
(Navarro et al, 1996; Navarro et al, 1997; Burkert, 1995). 

The gas component was assumed to be initially cold, Tq = 10^ K. As we 
see in our calculations, the influence of chaotic motions essentially reduces 
the dependence of model parameters on the adopted temperature cutoff 
and, therefore, on the adopted form of cooling function itself. 

The gas was assumed to be involved into the Hubble flow {Hq = 65 
km/s/Mpc) and into the solid - body rotation around z - axis. We added 
the small random components of velocities (A | v |= 10 km/s) to account 
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for the chaotic motions of fragments. The initial velocity field was defined 
as: 



v(a;, y, z) = [n{x, y, z) x r] + Hq ■ r + Av{x, y, z), (46) 

where il{x,y,z) is the angular velocity of an initially rigidly rotating 
system. 

The spin parameter in our numerical simulations is A pa 0.08. This 
parameter is defined in Peebles (1969) as: 



(47) 



G-iMgas + Mhalof/'^' 

Lq is the total initial angular momentum and Eq^ is the total initial 
gravitational energy of a protogalaxy. It is to be noted that for a system, 
in which angular momentum is acquired through the tidal torque of the 
surrounding matter, the standard spin parameter does not exceed A « 
0.11 (Steinmetz & Bartelmann, 1995). Moreover, its typical values range 
between A ^ 0.071^;^^, e.g. 0.02 < A < 0.11. 



3.2. DYNAMICAL MODEL 

At the final time step {tevol ~ 13.0 Gyr) the "star" and "gas" particle 
distributions have dimensions typical for a disk galaxy. The radial extension 
is approximately 25 — 30 kpc. The disk height is around 1 — 2 kpc. The "gas" 
particles arc located mainly within the central 5 — 10 kpc. 

We present in Fig. 4. the density distribution of gas p{r). Except for the 
central region (< 2 kpc), the gas distribution has a exponential form with 
radial scale length « 2.8 kpc. 

The column density distributions of gas (Tgas{r) and stars cr*(r) are 
presented in Fig. 5. The total column density is defined as: atot{^) = 
Cgas{f) + The total column density distribution atotir) is well ap- 

proximated (in the interval from 5 kpc to 15 kpc) with an exponential 
profile characterised with a « 3.5 kpc radial scale length. 

In the literature we found a large scatter of values for this radial scale 
length obtained for our Galaxy: 

- 5.5 kpc (van der Kruit, 1986), 

- 3.0 kpc (Eaton et al, 1984), 

- 2.5 kpc (Lewis &; Freeman, 1989), 

- 2.0 kpc (Jones et al, 1981), 

- 1.0 kpc (Mikami & Ishida, 1981), 

- 2.1 ± 0.3 kpc (Porcel et al, 1998), 

- 2.3 ± 0.1 kpc (Ruphy et al, 1996), 
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- 2.5toi kpc (Fux & Martinet, 1994). 

Our value is very close with a recently reported value: 3.5 kpc (Mera et 
al, 1998b). 

The value of atot ~ 55 Mq j-pc? near the location of the Sun (r « 9 kpc) 
is very close to a recent determination of the total density 52 it 13 M^/pc? 
(Mera et ol, 1998a). 

We present in Fig. 6. the rotational velocity distribution of gas Vrot{i'), 
at the end of calculation. In the figure we show the rotational curve for our 
Galaxy (Vallee, 1994). As we can see from the figure, the modelled disk 
galaxy rotation velocities are very close to our Galaxy rotation curve. 

The gaseous radial Vrad{f) and normal Vz{r) velocity distributions we 
present in Fig. 7. and Fig. 8. The radial velocity dispersion has a maximum 
value ~ 60 km/s in the center (this high vahic is mainly caused by the 
central strong bar structure). Near the Sun this dispersion drops down to 
K. 20 km/s. Such radial dispersion is reported in the kinematics study of the 
stellar motions in the solar neighborhood (Bienayme, 1998). The normal 
dispersion is near ~ 20 km/s in the whole disk. This value also coincides 
with the vertical dispersion velocity near the Sun (Bienayme, 1998). 

We present in Fig. 9. the temperature distribution of gas T{r). As we 
can see from the figure, the distribution of T(r) has a very large scatter from 
10^ K to lO*' K. Because, in our calculation we set the cutoff temperature 
for the cooling function at 10^ K, the gas can't cool to lower temperatures. 
The most (~ 90 %) of the gaseous mass in the galaxy have a temperatures 
~ 10^ K. The modeled process of SNII explosions injects to the gas a great 
amount of thermal energy and generates a very large temperature scatter. 
Such a temperature scatter is typical for our Galaxy's ISM. 

Therefore, as we can see, even our crude numerical approximation gives 
a good fit for all dynamical and thermal distributions of gas and stars in a 
typical disk galaxy like our Galaxy. 

3.3. CHEMICAL CHARACTERISTICS 

We present in Fig. 10. the time evolution of the SFR in galaxy SFR{t) = 
dM^{t)/dt. About 90% of gas is converted into stars at the end of calcu- 
lation. The most intensive SF burst happened in the first ~ 1 Gyr, with a 
maximal SFR f» 35 Mq/ut. After 1.5 — 2 Gyr the SFR is coming down 
like an "exponential function". The SFR has a value « 1 Mq/ut at the end 
of the simulation. 

To check the SF and chemical enrichment algorithm in our SPH code, 
we use the chemical characteristics of the disk in the "solar" cylinder (8 
kpc < r < 10 kpc). 




Figure 5. crgas{r), (Tt,{r) and atot{r) = (ygas{r)+<T*{r). The column density distribution 
in the disk of gas and stars in the final step 
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Figure 6. K-ot(f). The rotational velocity distribution of gas in the final step 
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Figure 7. Vradix). The radial velocity distribution of g£is in the final step 
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Figure 8. 14 (r). The perpendicular to disk normal velocity distribution of gas in the 

final step 
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Figure 9. T{r). The temperature distribution of gas in the final step 
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The age - metalicity relation of the " star" particles in the " solar" cylin- 
der, [Fe/H](t), we present in Fig. 11. The observational data in this figure 
are taken from Meusinger et al. (1991) and Edvardsson et al. (1993). 

We presented in Fig. 12. the metalicity distribution of the "star" parti- 
cles in the "solar" cylinder A^,,([Fe/H]). The model data are scaled to the 
observed number of stars (Edvardsson et al, 1993). 

The [0/Fe] vs. [Fe/H] distribution of the "star" particles in the "solar" 
cylinder we presented in Fig. 13. In this figure we also present the observa- 
tional data from Edvardsson et al. (1993) and Tomkin et al. (1992). 

All these model distributions are in good agreement, not only with pre- 
sented observational data, but also with other data collected from Portinari 
et al. (1997). 

The [0/H] radial distribution [0/H](r) we present in Fig. 14. The ap- 
proximation presented in the figure is obtained by a least - squares linear 
fit. At distances 5 kpc < r < 11 kpc the model radial abundance gra- 
dient is —0.06 dex/kpc. In the literature we found different values of this 
gradient defined in objects of different types. From observations of HII 
regions (Pcimbcrt, 1979; Shaver et al, 1983) wc obtained that the oxygen 
radial gradient is —0.07 dex/kpc. From observations of PN of different types 
(Maciel & Koppen, 1994) we obtained the values: —0.03 dex/kpc for PNI, 
-0.069 dex/kpc for PNII, -0.058 dex/kpc for PNIII, -0.062 dex/kpc for 
PNIIa, —0.057 dex/kpc for PNIIb. Therefore our model agrees well with 
the oxygen radial gradient in the our Galaxy. 

3.4. CONCLUSION 

This simple model provides a reasonable, self - consistent picture of the pro- 
cess of galaxy formation and star formation in the galaxy. Our calculations 
show that even a relatively small number of " gas" particles with physically 
motivated star formation criteria can reproduce the most of observed dy- 
namical and chemical characteristic in the galactic disks. The dynamical 
and chemical evolution of the modeled disk - like galaxy is coincident with 
the results of observations for our own Galaxy. The results of our modeling 
give a good base for a wide use of the proposed SF and chemical enrichment 
algorithm in other SPH simulations. 
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Figure 10. SFR{t) = dM*{t)/dt. The time evolution of the SFR in galaxy 
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Figure 11. [Fe/H](t). The age metalicity relation of the "star" particles in the "solar" 
cylinder (8 kpc < r < 10 kpc) 
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Figure 12. A^, ([Fe/H]). The metalicity distribution of the "star" particles in the "solar" 
cylinder (8 kpc < r < 10 kpc) 
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Figure 13. The [O/Fe] vs. [Fe/H] distribution of the "star" particles in the "solar" 
cylinder (8 kpc < r < 10 kpc) 
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